% ***************************************************************************************
% m_calibration_robust_popdens36.m 
% ***************************************************************************************
clear all; 
close all;
format shortG;

%maxNumCompThreads(1); % Uncomment for determinism stress-test; not required for reproducibility
% ==================== SET FUNDAMENTAL PARAMETERS ========================
gravity = 8;             % gravity coefficient (estimation in step2 gravity)   
rho = 0.66;              % discount factor (after 1955) 
rho50 = 0.66;            % --- (for 1950) 
sigma = rho/gravity;     % shape parameter of idiosyncratic shocks (after 1955)  
sigma50 = rho50/gravity; % --- (for 1950)

mu=0.15;                 % expenditure share on floor spaces (consumers) 
gammaL=0.7;              % cost share on labor (production)
gammaH=0.1;              % cost share on floor spaces (production)
eta=4;                   % floor space supply elasticity

% grid of alpha, beta 
lb = -0.25; 
ub = 0.50; 
gridint = 0.001;  
alphagrid=(lb:gridint:ub);  
betagrid=(lb:gridint:ub);      
Ngrid=size(alphagrid,2); 

% ========================= LOAD DATA  =========================================  
datafile = '../../data/processed data/quant/Hiroshima_DATA.csv'; CITYDATA = readtable(datafile);

ID = CITYDATA.geocode;   IDunique = unique(ID); 
if size(ID,1) ~= size(IDunique,1) 
    disp(' >>> CHECK DUPLICATED LOCATION DATA <<<');  pause 
end 
clear IDunique 
N=size(ID,1); Area = CITYDATA.area; GEOGRAPHY=[ID Area]; cbddist=CITYDATA.cbddist; 

Rdata36=CITYDATA.pop36; Rdata45=CITYDATA.pop45bomb; 
Rdata51=CITYDATA.pop51; Rdata55=CITYDATA.pop55; 
Rdata60=CITYDATA.pop60; Rdata65=CITYDATA.pop65;     
Rdata70=CITYDATA.pop70; Rdata75=CITYDATA.pop75; 
Rdata50=CITYDATA.pop50; Rdata53=CITYDATA.pop53;
Rdata57 = 0.6*Rdata55+0.4*Rdata60; 
Rdata66 = 0.8*Rdata65+0.2*Rdata70;

Ldata38=CITYDATA.emp38; Ldata45=CITYDATA.emp45bomb; 
Ldata53=CITYDATA.emp53; Ldata57=CITYDATA.emp57; 
Ldata63=CITYDATA.emp63; Ldata66=CITYDATA.emp66;     
Ldata75=CITYDATA.emp75; 

% Read latitude and longitude (for readability of the csv output in GIS)
lat = CITYDATA.lat; lon = CITYDATA.lon;

Ldata50=Ldata53*(sum(Rdata50)/sum(Rdata53)); 
Ldata55=0.5*Ldata53+0.5*Ldata57;  
Ldata60=0.5.*Ldata57+0.5.*Ldata63; 
Ldata65=(1/3)*Ldata63 + (2/3)*Ldata66; 
Ldata70=(5/9).*Ldata66+(4/9).*Ldata75; 

Rdata45=(sum(Ldata45)/sum(Rdata45)).*Rdata45;  
Rdata50=(sum(Ldata50)/sum(Rdata50)).*Rdata50;
Rdata55=(sum(Ldata55)/sum(Rdata55)).*Rdata55;  
Rdata60=(sum(Ldata60)/sum(Rdata60)).*Rdata60;
Rdata65=(sum(Ldata65)/sum(Rdata65)).*Rdata65;  
Rdata70=(sum(Ldata70)/sum(Rdata70)).*Rdata70;
Rdata75=(sum(Ldata75)/sum(Rdata75)).*Rdata75; 

% data for the post-war economy (every 5 years from 1950);
% NOTE: we use data from 1950 to 1975 for estimation
Rdata=[Rdata50 Rdata55 Rdata60 Rdata65 Rdata70 Rdata75];
Ldata=[Ldata50 Ldata55 Ldata60 Ldata65 Ldata70 Ldata75]; 
totalR = sum(Rdata,1); totalL = sum(Ldata,1);  
Popdata=totalR; 
if min(totalR ~= totalL) == 1   % check total pop ; 
    disp('>>> CHECK TOTAL POPULATION AND EMPLOYMENT ! <<< ');   pause
end 

% all data (every 5 years, 1945-75) 
Rdata_all = [Rdata45 Rdata50 Rdata55 Rdata60 Rdata65 Rdata70 Rdata75];
Ldata_all = [Ldata45 Ldata50 Ldata55 Ldata60 Ldata65 Ldata70 Ldata75]; 

% population/employment density; 
NT=size(Rdata,2);   
Rdensdata = Rdata./repmat(Area,1,NT);  
Ldensdata=Ldata./repmat(Area,1,NT);  

% COMPUTE LOWER BOUNDS for Theta; 
Rthetalb = (Rdata(:,2:NT)-Rdata(:,1:NT-1) < 0);  
Rthetalb = (1 - Rdata(:,2:NT)./Rdata(:,1:NT-1)).*Rthetalb;  
Lthetalb = (Ldata(:,2:NT)-Ldata(:,1:NT-1) < 0);  
Lthetalb = (1 - Ldata(:,2:NT)./Ldata(:,1:NT-1)).*Lthetalb; 
Rthetalb = max(Rthetalb); Lthetalb = max(Lthetalb); 
thetalb = ceil(100.*max([Rthetalb; Lthetalb]))./100;

% ================= LOAD FLOORSPACE PRICES ================================ 
l_q50=CITYDATA.lnQ_50; q50=exp(l_q50*(4/eta)); 
l_q55=CITYDATA.lnQ_55; q55=exp(l_q55*(4/eta)); 
l_q60=CITYDATA.lnQ_60; q60=exp(l_q60*(4/eta)); 
l_q65=CITYDATA.lnQ_65; q65=exp(l_q65*(4/eta)); 
l_q70=CITYDATA.lnQ_70; q70=exp(l_q70*(4/eta));
l_q75=CITYDATA.lnQ_75; q75=exp(l_q75*(4/eta)); 
Qdata=[q50 q55 q60 q65 q70 q75]; % Floor space price matrix  

% ================= LOAD COMMUTING COST ===================================
load('commuting_cost.mat', 'Kappa_mat');  

% ================= SET Theta VECTOR ======================================
thetavec=[0.53; 0.53; 0.53; 0.53; 0.53];

if min(thetavec - thetalb')<0
    disp('>>> CHECK THETA VALUES (LOWER BOUND) ! <<< ');   pause
end 

%  =============  GROUP TRACTS FOR GMM BY DISTANCE FROM CBD ===============
%  =============  N=174 BLOCKS to 5 GROUPS (BINS) ==>> 35 x 4 + 34 ======== 
% INDEX LOCATIONS BY ASCENDING ORDER OF CBD distance (FROM SMALL TO LARGE)
% sortn = location index when sorting locations by density from small to large
G=5;
Rpopdens36 = log(Rdata36./Area); 
[Gmat] = gmmiv(Rpopdens36,G,N);  GG = sum(Gmat,2); 

%  ========================================================================
%  ==================  MAIN INVERSION OF THE MODEL   ======================
%  ========================================================================
options0 = optimset('Display','none','Algorithm','trust-region-dogleg','MaxFunEvals',50000000, 'MaxIter',1000000,'TolFun',1e-6,'TolX',1e-6);
Xmat=zeros(N,NT); 
Ymat=zeros(N,NT); 
Kmat = zeros(N,N,NT); 
LQmat=zeros(N,NT); 

%  construct kernel for previous periods NT-1, ..., 1;  
Kmat(:,:,NT)=Kappa_mat(:,:,NT).^(-rho/sigma);   % last period 
for t=1:NT-1
    tt=NT-t; 
    theta=thetavec(tt); 
    if t < NT-1; Kmat(:,:,tt)=(Kappa_mat(:,:,tt).^(-rho/sigma)).*(Kmat(:,:,tt+1).^(rho*(1-theta))); 
    else;        Kmat(:,:,tt)=(Kappa_mat(:,:,tt).^(-rho50/sigma50)).*(Kmat(:,:,tt+1).^(sigma*rho50*(1-theta)/sigma50)); end
end 

% CONSTRUCT CONTINUATION VALUES FOR FLOOR SPACE PRICES 
LQmat(:,NT)=Qdata(:,NT);
for t=1:NT-1 
    tt=NT-t;
    theta=thetavec(tt); 
    if t < NT-1; LQmat(:,tt)=Qdata(:,tt).*(LQmat(:,tt+1).^(rho*(1-theta))); 
    else;        LQmat(:,tt)=Qdata(:,tt).*(LQmat(:,tt+1).^(rho50*(1-theta))); end
end 

K50=Kmat(:,:,1);
Q50=LQmat(:,1); 

%  =================== INVERSION OF FUNDAMENTALS =========================== 
for t=1:NT-1
    tt = NT-t;    % compute from the last period (1975,1970,1965,1960,1955)
    theta=thetavec(tt); 
    K = Kmat(:,:,tt+1); 
    Q=LQmat(:,tt+1);
    Rpre=Rdata(:,tt); Rpost=Rdata(:,tt+1); 
    Lpre=Ldata(:,tt); Lpost=Ldata(:,tt+1);
       
    X0=ones(N,1);  Y0=ones(N,1);     
    Rsyst = @(x)fReval2(x,K,Q,Rpre,Rpost,Lpre,Lpost,rho,sigma,theta,mu,gammaL,gammaH);
    Lsyst = @(y)fLeval2(y,K,Q,Rpre,Rpost,Lpre,Lpost,rho,sigma,theta,mu,gammaL,gammaH);  
    [x_fsolve, xfval_fsolve]=fsolve(Rsyst, X0, options0); X=abs(x_fsolve); X=X./geomean(X);  
    [y_fsolve, yfval_fsolve]=fsolve(Lsyst, Y0, options0); Y=abs(y_fsolve); Y=Y./geomean(Y);
    Xmat(:,tt+1)=X;    Ymat(:,tt+1)=Y;
end 
X55 = Xmat(:,2); Y55 = Ymat(:,2); 
X60 = Xmat(:,3); Y60 = Ymat(:,3); 
X65 = Xmat(:,4); Y65 = Ymat(:,4); 
X70 = Xmat(:,5); Y70 = Ymat(:,5); 
X75 = Xmat(:,6); Y75 = Ymat(:,6);

%  ===================== DECOMPOSE FUNDAMENTALS ===========================
amat3D=zeros(N,NT-1,Ngrid); 
bmat3D=zeros(N,NT-1,Ngrid);
a_errmat3D = zeros(N,NT-1,Ngrid); 
b_errmat3D=zeros(N,NT-1,Ngrid);
for ii=1:Ngrid
  
    alpha = alphagrid(ii); beta = betagrid(ii);
    [a,a_err] = ffund2(alpha,Ymat,Ldensdata,rho,sigma,thetavec,N,NT);     
    [b,b_err] = ffund2(beta,Xmat,Rdensdata,rho,sigma,thetavec,N,NT);

    amat3D(:,:,ii)=a;  a_errmat3D(:,:,ii)=a_err; 
    bmat3D(:,:,ii)=b;  b_errmat3D(:,:,ii)=b_err; 

end

%  ===========  GRID SEARCH: 1st Step GMM OBJECTIVE FUNCTION ===============
objvala = zeros(Ngrid,1);  
objvalb = zeros(Ngrid,1);
d_a_err = a_errmat3D(:,2:NT-1,:)-a_errmat3D(:,1:NT-2,:);   
d_b_err = b_errmat3D(:,2:NT-1,:)-b_errmat3D(:,1:NT-2,:);

for ii=1:Ngrid
    objvala(ii)=fobjm(Gmat,GG,d_a_err(:,:,ii),N,NT,G,eye(G*(NT-2))); 
    objvalb(ii)=fobjm(Gmat,GG,d_b_err(:,:,ii),N,NT,G,eye(G*(NT-2)));
end

[~, ii_mina]  = min(objvala);
[~, ii_minb]  = min(objvalb);

alphaest1st = alphagrid(ii_mina);
betaest1st = betagrid(ii_minb); 
[alphaest1st, betaest1st]

% ======================================================================
% ==========================  2nd Step GMM   ===========================
% ======================================================================
% Change in log residual fundamentals from first step, N x (NT-2) ; 
d_a_err1st = d_a_err(:,:,ii_mina);
d_b_err1st = d_b_err(:,:,ii_minb); 

mmadist1st=zeros(G,N,NT-2); mmbdist1st=zeros(G,N,NT-2);
for t=1:NT-2
    mmadist1st(:,:,t)=Gmat.*repmat(d_a_err1st(:,t)',G,1);  mmbdist1st(:,:,t)=Gmat.*repmat(d_b_err1st(:,t)',G,1);   
end
va = zeros(G*(NT-2),G*(NT-2)); vb=va; 
for i=1:N
    mmadistv = reshape(squeeze(mmadist1st(:,i,:)), G*(NT-2),1);  
    mmbdistv = reshape(squeeze(mmbdist1st(:,i,:)), G*(NT-2),1);  
    va = va + mmadistv*mmadistv'./N;   vb = vb + mmbdistv*mmbdistv'./N; 
end

objvala2=zeros(Ngrid,1);  
objvalb2=zeros(Ngrid,1);
wa = inv(va); wb=inv(vb);    % weight matrix 
for ii=1:Ngrid
    objvala2(ii)=fobjm(Gmat,GG,d_a_err(:,:,ii),N,NT,G,wa);
    objvalb2(ii)=fobjm(Gmat,GG,d_b_err(:,:,ii),N,NT,G,wb); 
end

[~, ii_mina2] = min(objvala2);
[~, ii_minb2] = min(objvalb2);

alphaest = alphagrid(ii_mina2);
betaest = betagrid(ii_minb2); 
[alphaest betaest] 


% ----- choose explicit, scale-aware finite-difference steps -----
h_a = max(1e-6, sqrt(eps)) * max(1, abs(alphaest));
h_b = max(1e-6, sqrt(eps)) * max(1, abs(betaest));

alphaest_se = alphaest + h_a;   
betaest_se  = betaest  + h_b;    

% compute moment function for parameter perturbation: 
[~,a_err_se] = ffund2(alphaest_se,Ymat,Ldensdata,rho,sigma,thetavec,N,NT);  % N x (NT-1): period 2,3,...,T
[~,b_err_se] = ffund2(betaest_se,Xmat,Rdensdata,rho,sigma,thetavec,N,NT);    
d_a_err_se=a_err_se(:,2:NT-1)-a_err_se(:,1:NT-2);    
d_b_err_se=b_err_se(:,2:NT-1)-b_err_se(:,1:NT-2);     

mma_se=zeros(G,N,NT-2); 
mmb_se=zeros(G,N,NT-2); % G x N x NT-2
for t=1:NT-2
   mma_se(:,:,t) = Gmat.*repmat(d_a_err_se(:,t)',G,1);  
   mmb_se(:,:,t) = Gmat.*repmat(d_b_err_se(:,t)',G,1);  
end

[aa,a_err] = ffund2(alphaest,Ymat,Ldensdata,rho,sigma,thetavec,N,NT) ; 
[bb,b_err] = ffund2(betaest,Xmat,Rdensdata,rho,sigma,thetavec,N,NT) ;   
d_a_err = a_err(:,2:NT-1)-a_err(:,1:NT-2);    
d_b_err = b_err(:,2:NT-1)-b_err(:,1:NT-2);    
mma=zeros(G,N,NT-2); mmb=zeros(G,N,NT-2);% G x N x NT-2
for t=1:NT-2
    mma(:,:,t)=Gmat.*repmat(d_a_err(:,t)',G,1); 
    mmb(:,:,t)=Gmat.*repmat(d_b_err(:,t)',G,1);   % G x N
end

% ----- build numerical gradients (stacked moments), divide by actual steps -----
grad_a = squeeze(sum(mma_se - mma, 2));   % G x (NT-2)
grad_b = squeeze(sum(mmb_se - mmb, 2));   % G x (NT-2)

% Stack to vectors and normalize by N 
grad_a = reshape(grad_a, G*(NT-2), 1) / h_a / N;   % M x 1, M = G*(NT-2)
grad_b = reshape(grad_b, G*(NT-2), 1) / h_b / N;

% ----- make weight matrices symmetric (safer numerics) -----
va = (va + va.')/2;
vb = (vb + vb.')/2;

% ----- compute scalar GMM variances without inv() -----
s_a = grad_a' * (va \ grad_a);     
s_b = grad_b' * (vb \ grad_b);     
v_alpha_gmm = sqrt(1 / (s_a * N));
v_beta_gmm  = sqrt(1 / (s_b * N));

% =================== SAVE PARAMETER RESULTS ========================
objresult=[mu,gammaL,gammaH,eta,gravity,rho,rho50,sigma,sigma50,alphaest,betaest,v_alpha_gmm,v_beta_gmm];
resultparam=array2table(objresult,'VariableNames',{'mu','gammaL','gammaH','eta','gravity','rho','rho50','sigma', 'sigma50', ... 
    'alpha','beta','sdalpha','sdbeta'});  
writetable(resultparam,'../../output/quant/table_parameter_estimates_popdens36.csv'); 

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%                       FUNCTIONS
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ---------- FUNCTION FOR CONSTRUCTING IV --------- 
function[gmat]=gmmiv(ivvar,nG,N); 
% ivvar: variable for constructing IV
Ground=round(N/nG); gn=zeros(N,1); gmat=zeros(nG,N);
[~,sortn] = sort(ivvar,'ascend'); 
Gthres=[Ground+1; 2*Ground+1; 3*Ground+1; 4*Ground+1; 5*Ground+1; 6*Ground+1; 7*Ground+1; 8*Ground+1; 9*Ground+1];
for i=1:N
    ii=find(sortn==i);
    if ii<Gthres(1); gn(i)=1; 
    elseif ii >= Gthres(nG-1); gn(i)=nG; 
    else; 
        for g=2:nG-1 
        if Gthres(g-1) <= ii && ii < Gthres(g); gn(i)=g;  
        else; continue  
        end 
        end 
    end 
end
for i=1:N
    gi = gn(i);
    gmat(gi,i)=1;
end
end 

% --------  FUNCTION FOR COMPUTING CHARACTERISTICS BASED ON GRAVITY STURCTURE --------
function[Reval] = fReval2(X,K,Q,DRpre,DRpost,DLpre,DLpost,rho,sigma,theta,mu,gammaL,gammaH)   % eq D.16. X corresponds to Theta

N = size(DRpre,1); 
X=abs(X);  % avoid negative values
X=X./geomean(X);
%xtemp= K.*(repmat(X',N,1).^(rho/sigma)); 
xtemp=K.*((repmat(Q',N,1).^(-mu)).*repmat(X',N,1)).^(rho/sigma);
xtemp=xtemp./repmat(sum(xtemp,2),1,N); 

Rdif=DRpost-(1-theta).*DRpre; 
Ldif=DLpost-(1-theta).*DLpre; 

Reval = Rdif - sum(xtemp.*repmat(Ldif,1,N),1)'; 
Reval = abs(Reval); 

end

function[Leval] = fLeval2(Y,K,Q,DRpre,DRpost,DLpre,DLpost,rho,sigma,theta,mu,gammaL,gammaH) % eq D.17. Y corersponds to Omega

N = size(DRpre,1); 
Y=abs(Y);  
Y=Y./geomean(Y); 
ytemp= K.*((repmat(Q,1,N).^(-gammaH/gammaL)).*(repmat(Y,1,N).^(1/gammaL))).^(rho/sigma);
ytemp=ytemp./repmat(sum(ytemp,1),N,1); 

Rdif=DRpost-(1-theta).*DRpre; 
Ldif=DLpost-(1-theta).*DLpre; 

Leval = Ldif - sum(ytemp.*repmat(Rdif',N,1),2); 
Leval = abs(Leval); 

end

% --------  FUNCTION FOR COMPUTING THE EXOGENOUS FUNDAMENTALS --------
function[c, c_err] = ffund2(param,Zmat,densdata,rho,sigma,thetavec,N,NT)  

Z0=Zmat(:,NT); cmat=zeros(N,NT-1);
for t=1:NT-1
    denspost=densdata(:,NT-t+1); 
    invspill = denspost.^(-param);
    cc = invspill.*Z0;
    cmat(:,NT-t)= cc; 
    theta=thetavec(NT-t); 
    Z0=Zmat(:,NT-t).*(Z0.^(-rho*(1-theta)));     
end

c = log(cmat);   % N x NT-1
c_te = sum(c,1)./N;
c_err = c-repmat(c_te,N,1); 
end 

% ----------  FUNCTION FOR OBJECTIVE FUNCTION ------------
% Gmat: G (number of grids) x N (number of locations)  
% cmat: N (number of locations) x NT-2 (number of periods well defined for changes)
% Wmat: weight matrix, G*NT-2 (number of grids) x G*NT-2 
% objdist: scalar 
function[objdist] = fobjm(Gmat,GG,cmat,N,NT,G,Wmat); 
mdist = zeros(G,N,NT-2);
for t=1:NT-2
    mmt = Gmat.*repmat(cmat(:,t)',G,1);  % G x N
    mdist(:,:,t) = mmt;                 
end
mdist = squeeze(sum(mdist,2));   %  G x NT-2
mdist = reshape(mdist, G*(NT-2),1); 
objdist=(mdist./N)'*Wmat*(mdist./N);
end
